/*
This master file should be able to 
1) prepare the data for the 36 cell case
2) fit various models to this data (specifically JD, honesty, JD+honesty, DD and KS)
3) calculate the AIC
4) make a nice table to summarise the results of this work. 


Note for ks and dd this calls other do files, which fit these models in the 
*/

*****************************************************************
* 1) prepare the data for the 36 cell case
*****************************************************************

* data prep, for the 36 cells
clear all

forvalues i = 21 / 22 {
use "data/both.dta", clear
keep if treat==`i'
drop treatment treat experiment
gen reports_`i'=1
collapse (sum) reports_`i', by(roll1 roll2) 
save  "temp/t`i'.dta", replace
}

clear
set obs 36
// Generate roll combinations (0,0), (0,1), ..., (5,5)
gen roll1 = floor((_n-1)/6)
gen roll2 = mod(_n-1, 6)

joinby using "temp/t21.dta", unmatched(both)  _merge(m21) 
replace reports_21=0 if reports_21==.

joinby using "temp/t22.dta", unmatched(both)  _merge(m22)
replace reports_22=0 if reports_22==.


* turn into fractions
sum reports_21
di `r(sum)' 
replace reports_21=reports_21/`r(sum)' 

sum reports_22
di `r(sum)' 
replace reports_22=reports_22/`r(sum)' 

*sum 
sum 
drop m21 m22

save "temp/36cell.dta", replace

use "temp/36cell.dta", clear


*****************************************************************
* 2) fit various models to this data (specifically JD, honesty, JD+honesty, DD and KS)
** let's start with JD, treatment 22
*****************************************************************


* Create a new frame to store results
frame change default
capture frame drop results
frame create results
frame change results

* Set up the structure
clear
set obs 0
gen treatment = .
gen str10 model = ""
gen aic = .
gen bic = .
gen rmse = .

* back to current frame
frame change default

* now use a local to select treatment. 
gen reports=.
gen predicted_honesty=1/36
gen predicted_jd=(roll1==roll2)/36
replace predicted_jd=2/36 if roll1>roll2

forvalues j=21/22 {
global treat `j'
replace reports=reports_$treat

global model_name "KS"
do "code/ks.do"

global model_name "DD"
do "code/dd.do"
}

foreach var of varlist predicted_* reports_* {
 replace `var' = `var'*100
 }

// Clear any existing stored estimates
estimates clear

forvalues i = 21/22 {
    foreach var of varlist predicted_honesty predicted_jd predicted_dd_`i' predicted_ks_`i' {
           // Set parameter number for IC calculation
        if inlist("`var'", "predicted_honesty", "predicted_jd") {
            local parameter_number 0
        }
        else if inlist("`var'", "predicted_dd_`i'") {
            local parameter_number 1
        }
        else {  // predicted_ks_`i'
            local parameter_number 2
        }
		
		 // make name readable
    local clean_var `=subinstr("`var'", "predicted_", "", 1)'_`i'

        constraint 1 `var'=1
        qui: cnsreg reports_`i' `var' , constraint(1) nocons
        scalar est_rmse  =  e(rmse) * sqrt(e(df_r) / e(N))

		estat ic, df(`parameter_number') n(36) 
		matrix ic_results = r(S) 
		scalar aic_value = ic_results[1,5]
		scalar bic_value = ic_results[1,6]
		* Store in results frame
		frame change results
		set obs `=_N+1'
		replace treatment = `i' if _n==_N
		replace model = "`clean_var'" if _n==_N
		replace aic = aic_value  if _n==_N
		replace bic = bic_value  if _n==_N
		replace rmse =  est_rmse if _n==_N
		frame change default

    }
}

frame change results
order model treatment aic bic rmse
replace aic=round(aic,.01)
replace bic=round(bic,.01)
replace rmse=round(rmse,.001)
list
xpose, clear varname
order _
drop if _n<.5

rename v1 Honest_2s
rename v2 JD_2s
rename v3 DD_2s
rename v4 KS_2s
rename v5 Honest
rename v6 JD
rename v7 DD 
rename v8 KS

drop if _n==1
list
texsave using "output/table4.tex", replace  frag



